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Abstract. The mechanism of cold- and pressure-denaturation are matter of debate. 
Some models propose that when denaturation occurs more hydrogen bonds between 
the molecules of hydration water are formed. Other models identify the cause in the 
density fluctuations of surface water, or the destabilization of hydrophobic contacts 
because of the displacement of water molecules inside the protein, as proposed for high 
pressures. However, it is clear that water plays a fundamental role in the process. Here, 
we review some models that have been proposed to give insight into this problem. Next 
we describe a coarse-grained model of a water monolayer that successfully reproduces 
the complex thermodynamics of water and compares well with experiments on proteins 
at low hydration level. We introduce its extension for a homopolymer in contact with 
the water monolayer and study it by Monte Carlo simulations. Our goal is to perform 
a step in the direction of understanding how the interplay of cooperativity of water 
and interfacial hydrogen bonds affects the protein stability and the unfolding. 
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1. Introduction 

One of the most intriguing challenge in biological physics is the nature of protein folding- 
unfolding processes. The temperature range of stability of proteins is in general small. 
For example, staphylococcal nuclease (Snanse-a small protein containing 149 amino- 
acids) folds at low pressure is approximately between 260 K and 320 K pp. 

Heat destabilizes proteins. By increasing the bath temperature T, thermal 
fluctuations increase and disrupt the folded configurations of proteins. Usually, by 
decreasing T, proteins crystallize, but surprisingly some proteins unfold at sufficient 
low temperature instead of crystallizing [U El El IU El El E] ■ Cold denaturation seems 
to be a general phenomenon for proteins, generally occurring well below 0°C, the 
freezing point of water. In same cases, for example for Snase pQ, the cold denaturation 
cannot be directly observed, but experimental data can be extrapolated to predicted 
the lower temperature of stability for the protein. More generally, to make the cold 
denaturation observable destabilizing agents can be used. Interestingly, Pastore et al. [2] 
observed that Yeast frataxin under physiological conditions undergoes cold denaturation 
below 7°C and remains folded up to 30°C Hence, Yeast frataxin could be an excellent 
prototype for studying folding-unfolding transition, both hot and cold, under accessible 
conditions. 

Proteins can unfold also by pressurization. It has been observed that the increase 
of pressure induces the unfolding of protein [H [9]. The pressure-unfolding process can 
be rationalized by considering that the folded structure usually includes cavities. High 
pressure can induce elastic response of the protein, deforming its structure and pushing 
water molecules inside the cavities. The water molecules from inside would swell the 
protein, with consequent loss of protein functionality [8]. Because is difficult to to 
separate the protein response to high hydrostatic pressure from the response of the 
aqueous environment, the understanding of the problem is still under debate. 

1.1. Thermodynamics of proteins unfolding 

By increasing the thermal energy k^T [k^ is the Boltzmann constant), the protein 
residues vibrate faster, accessing new possible configurations, i. e. increasing the the 
entropy S of the system. This increase leads to the hot denaturation, in the same way 
an increase of ksT leads to the melting of a crystal. 

The cold denaturation instead, cannot be explained as the effect of an increase of 
entropy. By decreasing T, the entropy of the system decreases. This is why we cannot 
melt a crystal by cooling. Hence, in the case of proteins there must be a complex 
mechanism that induces the cold denaturation. To understand this mechanism is 
necessary to introduce the concept of Gibbs free energy G = H—TS, where H = U+PV 
is the the enthalpy of the system, U internal energy of the system, V the volume and 
P the pressure. 

General principles of thermodynamics tell us that at any value of T and P the 
system minimizes its Gibbs free energy, where the system, in our case, is the solution of 
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proteins and water. Hence, the free energy balance must take into account both water 
molecules and protein residues. The experimental fact that solvated proteins unfolds by 
decreasing T means that at lower T the difference 

AG = G U — Gf (1) 

between the unfolded (u) and folded (/) states is 

AG = AH - TAS < 0, (2) 

where AH = H u — Hf and AS = S u — Sf. 

The total variation of the entropy of the system is given by AS = AS P + AS W 
where AS P and AS W are the variation of entropy of the protein residues and water 
molecules, respectively. By unfolding, the protein entropy increases, AS P > 0. On the 
other hand, the protein contribution to AH is positive, AH P > 0, because the enthalpy 
of the protein increases when the protein unfolds (H p is proportional to the number of 
contact points of the protein). Therefore, the protein contribution to AG, AH p — TAS p , 
could be negative or positive depending on the relative variations and on T and does not 
guarantees that Eq. ^ is satisfied. Hence, water contribution to the total balance of 
Eq. ([2]) could be relevant. To date, is widely shared the idea that the native-folded state 
is stabilized by the quasi-ordered network of water molecules hydrating the non-polar 
monomers [?]. 



1.2. Protein Phase Diagram 

Experiments are consistent with a protein stability phase diagram with an elliptic shape 
EH El El SI El [10] in P-T plane (Fig. [TJ . Outside the elliptic region the protein unfolds 
loosing its biological function. Following Hawley [TH [12] , we can calculate AG of the 
whole system (protein and water) assuming that a protein can stay in only two distinct 
states, folded and unfolded as in Eq. ([I]). 

Because the internal energy of the system U is a state variable, we can express its 
infinitesimal variation dU as a function of two thermodynamical quantities. If we assume 
that the unfolding process can be described by infinitesimal quasi-static transformations, 
applying the First Law of Thermodynamics, we get 

dU = 6Q- 5W, (3) 

where 5Q and SW are the infinitesimal heat absorbed and the infinitesimal work done 
by system, respectively, along the generic transformation. Since at constant T and P is 
5Q = TdS and 5W = PdV, we can express the internal energy variation, for a constant 
number of particles N, as 

dU = TdS - PdV = dU{S,V). (4) 

Differentiating H = U + PV we get 

dH = dU + PdV + VdP = TdS + VdP = dH(S,P), (5) 
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Figure 1. Schematic representation of the phase diagram of a protein. Within 
the elliptic shape the protein is folded, while it unfolds by increasing temperature T 
(hot denaturation), by decreasing T (cold denaturation), by increasing or decreasing 
pressure P (pressure denaturation). Each folding-unfolding process is characterized 
by different variation of entropy AS and variation of volume AV. The axes of the 
ellipse are loci where AS = and AV = (see text for discussion). Adapted from [3]. 



and differentiating G = H — TS, we finally get 

dG = dH- TdS - SdT = -SdT + VdP = dG(T, P). 

Hence, it is 

dAG = -ASdT + AVdP 



(6) 



(7) 



with AS = S u — Sf and AV = V U — Vf. By expanding AS and AP to the first order 
around ASq and AVq, changes at T and Pq, we get 



A^AS 0+ (^) p( r-r o)+ (^) T (P-P„) 



(8) 
(9) 



and from Eq. ([7])-([9]), by integration, 



AG(P,T) = mP - P ) 2 + 2Aa(P - P )(T - T ) + 



(10) 



- AC P [(T - T ) - T ln(T/T )] + AV (P - P ) - AS (T - T ) + AG , 

where a = (dV/dT)p = —(dS/dP)T is the thermal expansivity factor, related to the 
isobaric thermal expansion coefficient apby ap = a/V; Cp = T(dS/dT) p is the isobaric 
heat capacity and (3 = (dV/dP)T is the isothermal compressibility factor related to the 
isothermal compressibility Kt by the relation Kt = —(/3/V). All the quantities with 
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the subscript equal to zero are usually referred to ambient conditions. By developing 
the logarithm to the second order around (T , Pq) 

T\ T - T (T — T ) 2 
To) ~ 27f 



(ID 



we get 



AG(P, T) = f(P- P ) 2 + 2Aa(P - P Q )(T - T )+ 



(12) 



- -27b" ( T " T °) + A ^( P " P °) " AS ^ T - + AG « 

that is the equation of en ellipse given the constraint 



Aa 2 > ACpAf3/T . (13) 

This condition is guaranteed by the different sign of ACp and A/3, as can be observed 
for some proteins, as reported by Hawley [TT] . 

The Eq. (12) is AG(P,T) Taylor expansion arrested to the second order, holding 
for Aa, A/3 and ACp independent of T and P, as generally valid. Adding third order 
terms in the expansion makes minimal effects on the elliptic shape of the stability region. 

At maximum pressure -P max of stability for the protein, dAG/dT = AS = 0, while 
at the maximum temperature T max of stability, dAG/dP = AV = 0. Therefore, based 
on Hawley's theory it is possible to make general predictions about the changes of AV 
and AS" as schematically summarized in Fig. [TJ This phenomenological theory has no 
explicit information about the protein structure, and makes strong assumptions, such as, 
for example, that the protein only has two states, or that equilibrium thermodynamics 
holds during the denaturation. The last assumption, in particular, implies that the 
all process would be reversible. Nevertheless, consistency with Hawley's theory is a 
good test for models of protein unfolding. In the next section we review some of these 
models. The review does not pretend to be exhaustive, but it has the aim of mentioning 
a number of positive results of the theory of protein folding. 



2. Models for protein unfolding 

In 1989 Lau and Dill proposed the HP model for protein folding [13J. By assuming 
that the exposed surface of hydrophobic residues is energetically unfavorable at low T, 
the model reproduces the folding of the protein (hydrophobic collapse). The protein is 
represented as a self- avoiding chain on a lattice. The chain is composed by two different 
categories of amino acids: H (non-polar) and P (polar). The presence of the aqueous 
environment is taken into account in an effective way, by introducing an attractive 
contact interaction between H monomers. No other interactions are present in the 
system. 

Under these hypothesis, the authors show that the features of the folding process 
depend on the HH energy interaction, the length of the chain, and the specific sequence 
of H and P monomers. Moreover, for long chains one folded state dominates. 
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The model has the virtue to reduce the complexity of the folding process to a 
manageable level. All the electrostatic and chemical properties of each amino-acid are 
simplified by allowing only two possible states. The degrees of freedom of the solvent 
are not explicitly included and the driving force for the folding is the hydrophobic 
interaction of non polar monomers. Nevertheless, the HP model cannot describe cold 
denaturation. Therefore, the experimental evidence of cold denatured proteins calls for 
a reconsideration of the hydrophobic interaction and its dependence on temperature and 
structure of hydration water [H El El Ej- 

Back in 1948, Frank and Evans [Hj discussed the tendency of water to form ordered 
structures around non polar solutes to minimize the free energy cost of solvation. 
As a consequence, hydrophobic solutes are "structure makers" for water, facilitating 
the formation of cages around the solute. The effect of these structures around 
hydrophobic solutes is to reduce the entropy with respect to the bulk and to compensate, 
approximately, the enthalpy cost for the creation of a cavity to allocate the solute. 

As discussed by Muller in 1990 [15], the compensation of the enthalpy implies that 
water-water hydrogen bonds (HBs) at the interface with the hydrophobic solute are 
stronger than those in the bulk. This is consistent with the experimental observation 
that the excess molar heat capacity for a nonpolar solute at infinite dilution in water is 
positive. This quantity, defined as the difference of the partial molar heat capacity in 
solution with the heat capacity of the pure liquid solute, is far larger at 25 °C when the 
solvent is water than for any other solvent [151 US] ■ 

The statement that HBs are stronger at the hydrophobic interface has led to the 
misconception that water around a hydrophobic solute has an iceberg-like structure. 
Computer simulations [T71 Q2] , theoretical analysis [121 EDI EI] , and neutral scattering 
studies [22J are inconsistent with iceberg-like structures. Hence, the restructuring of 
water around a solvent seems not to play a relevant role in the hydrophobic effect. 
Nevertheless, Muller [15] showed that if hydration HBs are enthalpically stronger but 
fewer than in bulk, a model with two-states HBs can reproduce the sign reversal of the 
proton NMR chemical shift with T and the heat capacity change upon hydration. 

On the other hand, a common opinion [231 EI] is that the large free-energy 
change associated to the hydrophobic effect is due to the small size of the water 
molecules with respect to the solutes, and that the free-energy change associated to the 
network reorganization around hydrophobic particles is negligible due to compensation 
of enthalpy and entropy, although it may account for the large heat capacity change upon 
hydration. This observation apparently ruled out Muller model, where the enthalpy- 
entropy compensation upon hydration was not present. 

Nevertheless, Lee and Graziano in 1996 [25] showed that Muller model can be 
slightly modified to recover also the enthalpy-entropy compensation upon hydration. 
The Muller-Lee-Graziano model was further simplified by De los Rios and Caldarelli in 
2000 [261 E3 EH] in order to reduce the number of parameter. By further simplifying the 
the description of bulk water, they recovered hot and cold denaturation for a protein 
represented as a hydrophobic homopolymer. A development of this model has been 
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recently used to study the effective interaction between chaotropic agents and proteins 
[29]. 

The model by De los Rios and Caldarelli has been generalized by Bruscolini and 
Casetti (30J ]3Tj in 2000 by allowing each monomer of non-polar homopolymer to be 
in contact with a cluster of water molecules. Each cluster has an infinite number of 
possible states and only one state minimizes the free-energy cost of the interaction 
with the hydrophobic monomers. The model reproduces the trends of thermodynamics 
averages in accordance with experiments [32] and simulations [33], and predicts the 
warm and cold denaturation. These results are qualitatively similar to those of the 
Muller-Lee-Graziano model, further supporting the relevant role of the solvent in the 
folding-unfolding process. 

Cold denaturation and T-dependence of the hydrophobic effect were also observed 
by Dias et al. analyzing a nonpolar homopolymer in Mercedes-Benz (MB) water |34j. 
The MB model, originally introduced by Ben-Nairn [35], represents water molecules 
as disks in two-dimensions with three possible HBs (arms) as in a Mercedes-Benz 
logo. Water molecules interact via van der Walls potential and HB interactions. HB 
interaction is modeled with a Gaussian potential, favoring a fixed value for the water- 
water distance and aligned arms for facing molecules. Simulations show that the average 
HB energy is higher for shell water than for bulk water at high T, while is lower at lower 
T. Therefore, by cooling the solution, is energetically more convenient to increase the 
protein surface exposed to water, inducing protein unfolding. In this model, the water 
molecules forming a cage around the protein monomers are strongly H-bonded to each 
other. The highly ordered structure of the solvent around the monomers decreases the 
entropy of water, compensating the increase of the entropy associated to the protein 
unfolding. 

This model has been criticized [30] because it assumes, without proof, that the 
enthalpy gain dominates at low T, giving rise to free-energy gain upon unfolding 
of the protein. In particular, Yoshidome and Kinoshita [36J analyzed by integral 
equation theory the behavior of a nonpolar homopolymer composed by fused hard- 
spheres of different diameters immersed in smaller hard spheres, with permanent 
electrostatic multiple moments, representing the solvent [37J. The protein-water 
interaction is represented by a hard sphere potential and water-water interaction by 
a hard sphere potential and an electrostatic contribution given by the electrostatic 
multipole expansion. The author found that denaturation is characterized by large 
entropy loss and large enthalpy gain. However, these two contributions to the free 
energy almost completely cancel out and make no significant contribution to the free- 
energy change. They found that the driving mechanism for cold denaturation is the 
translational entropic-loss of water due to the large excluded volume of the hydrophobic 
particles. They observed that at low T water diffuses less, therefore the hydrophobic 
effect is weaker and the protein unfolds. 
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2.1. Pressure effects 

Pressure effects on protein denaturation have been also considered in microscopic 
theoretical models. For example, in 2003 Marques and coworkers [38] considered a 
hydrophobic homopolymer, represented as self-avoiding random walk, embedded in a 
water bath on a compressible lattice model in two-dimensions. Water-water interactions 
are represented by the Sastry et al. water model [39]. The polymer-water interaction 
is repulsive, being proportional to the density number of HBs and to the number of the 
missed native contact points of the protein. The model displays warm denaturation, 
pressure denaturation and cold denaturation at high pressure in agreement with stability 
diagram of some proteins [10] , although not with others [H] . A peculiarity of this model 
is that the effective repulsion between protein and solvent is coupled to the average 
number of HBs of bulk water. 

To remove this coupling to an average property of the bulk, in 2007 Patel and 
coworkers [12] proposed a model where water at the interface with protein has a 
restricted number of accessible orientations for the HBs compared to the bulk. Along 
with the entropic cost, the interfacial HBs also have an additional enthalpic bonus with 
respect to bulk water, following the ideas discussed by Muller, Lee and Graziano. The 
model displays a stability phase diagram with hot, cold and pressure denaturation. 
However, it does not reproduce all the expected features the schematic phase diagram 
of Fig. 0. In particular, the model does not reproduce the elliptic shape of the phase 
diagram and the low-P region with AS > and AV > for hot denaturation. These 
results were confirmed by extending the model to the case with heteropolymers. 

In the attempt to reproduce the elliptic phase diagram for protein stability, we 
propose here a model starting from the assumption that HBs at the interface are stronger 
compared to HBs in bulk water [13]. We'll proceed as follows: in section 3 we describe 
the model for nano-confined water, in section 4 we summarize recent results for the 
model, to clarify its water-like behavior. In the section 5 we propose a protein-water 
interaction mechanism displaying some preliminary results. 

3. Hydrophobic nanoconfinement for water 

We consider a monolayer of water nano-confined between hydrophobic plates. The 
interaction between water molecules and the surface is represented by a hard-core 
repulsion. The confinement is such to inhibit the formation of bulk water structures. 
For example, bulk water is known to preferentially form four HBs with four nearest 
neighbor molecules in an approximate tetrahedral structure at low temperature and 
pressure [H]. Hard confinement inhibits the formation of such bulk structure. For 
example, Kumar et al. [15] found, by molecular dynamics simulations with TIP5P- 
water confined between flat hydrophobic plates separated by 0.7 nm, an almost-flat 
monolayer of water molecules, each with four neighbors in an orientationally-disordered 
square lattice. 
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To define a tractable model we coarse-grain this structure of water under similar 
conditions [3S1 El SHI HH] . We divide the volume between the hydrophobic plates, and 
accessible to water, into N cells each containing one water molecule and each with 
volume V/N = r 2 h, where h is the separation between the flat planes, r > 2ro is the 
average distance between water molecules, with r equal to the van der Waals radius 
of a water molecule. The van der Waals attraction (due to dispersion forces) and the 
repulsive interactions (due to the Pauli exclusion principle) between water molecules are 
described by a Lennard- Jones interaction 



12 

r \ r 



(14) 



where r^ is the distance between molecules i and j and the sum is performed over all 
the molecules. 

To each cell we associate a variable n, = 0, 1. If the number density pi of the cell i 
is pir 2 h > 0.5 then rij = 1, otherwise rij = 0. 

To take into account the decrease of orientational entropy due to the formation 
of HBs, we introduce for each water molecule i four bonding indices o^-, one for each 
possible HB. Each variable <7jj can assume q different values, = l...q. We choose 
the parameter q by selecting 30° as the maximum deviation from linear bond (i.e. 
q = 180°/30° = 6). Hence, every molecule has q A = 1296 possible orientations. 

The covalent (directional) HB attraction component is expressed by the 
Hamiltonian term 

^HB = - J ^2 "i".! (S " , <S ", - 

<ij> 

where J > represent the covalent energy gained per HB, the sum is over nearest 
neighbors cells, and 5 a & = 1 if a = b, otherwise. 

The experiments show that the formation of a HB leads to an open structure that 
induces an increase of volume per molecule [4*H 150] . This effect is incorporated in the 
model by considering that the total volume of the system is 

V = V + iWffiB, (16) 

where Vq is the volume of the system without HBs, and t>HB is the increment due to the 
HB. 

The term =^hb quantifies the two-body component for HB interaction. On the other 
hand, the distribution of 0-0-0 angle shows a strong T-dependence [51] that suggest 
the presence of many-body component for HB interaction. We quantify this component 
by the Hamiltonian term 

^Coop = - Jar ^ ^2 ( 17 ) 

i (k,l)i 

where Jo- > is the characteristic energy of the cooperative component. The sum is 
performed over the six different pairs (k, l)i of arms of the molecule i. Therefore the 
total Enthalpy for the water is 

H = U + Jt nB + je Coop + PV = U-(J- Pv UB )N RB - JN a , (18) 
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where 



Nsb = ^2 ">".i (S « . (S » 



(19) 



<ij> 

is the total number of HB and 



(20) 



i {k,l)i 



is the total number of HBs optimizing the cooperative interaction. 
4. Results for water model 

We study our model by Monte Carlo (MC) simulations and mean-field calculations. MC 
simulations are performed in the NPT ensemble where the volume of the system is a 
stochastic variable. We consider periodic boundary conditions in the directions parallel 
to the confining surfaces. 

4-1. Liquid-gas phase transition and anomalies. 

Previous calculations have shown that the system displays a liquid-gas first-order 
phase transition ending in a critical point C at approximately kpTc/e = 1.9 ± 0.1 
and Pcva/t = 0.80 ± 0.05 (IS, H?] SHI E2J E31 EI] , in qualitative agreement with mean 
field results HHSSHS]. 

The model reproduces several anomalies of water. For example, the system presents 
density anomaly, i.e. the isobaric increase of density upon cooling, up to reach a 
temperature of maximum density (TMD). The system also displays diffusion anomalies 
[55], maxima of isothermal compressibility Kt, isobaric heat capacity Cp and the 
isobaric thermal expansion coefficient atp (1H1 HHJ ESI EH EH E3] related to the anomalous 
behavior of water in the supercooled region. 

4-2. Dynamical slowing down of water in supercooled region. 

The dynamical behavior of the model at low T has interesting features ESI EZ] • The 
dynamics of the HBs at constant P displays an increase of the correlation time when T 
is decreased. The increase is faster at higher T then at lower T and shows a crossover at 
the temperature when Cp has a maximum [561 EZ]- Results clarify that the crossover is 
due to a structural change in the HBs network [561 EZ] • The qualitative features of this 
crossover have been successfully compared to experimental results for confined water at 
increasing P [59]. In particular Franzese and de los Santos [52] have been showed that 
at high pressure (P ~ 2000 bar) the effect of HB is negligible due to the high enthalpic 
cost to form a HB and the correlation function decays as an exponential. At low P 
(P ~ 1 bar) the correlation is large also at long times and the system get stuck in a 
glassy state. The structural analysis shows that under these conditions the HB network 
develops gradually by decreasing T and traps the system in metastable configurations. 
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For intermediate value of pressure the correlation function C(t) is well described by a 
stretched exponential function 

C(t) = C e-W, (21) 

where C , r and (3 < 1 are fitting constant ((3 = 1 correspond to exponential decay). 
The quantity 1 — (3 is a measure of the heterogeneity in the system. As we approach to 
a characteristic value of pressure P c > , [3 reaches its minimum value ((3 ~ 0.4). This is 
consistent with the experimental value of (3 ~ 0.35, observed for intermediate scattering 
correlation function of water hydrating myoglobin at low hydration level (h = 0.35 g 
H 2 0/g protein) [501 EI]- Therefore, Franzese and de los Santos result indicates that the 
system exhibits a largest amount of heterogeneity at Pc>- As we will discuss in the next 
section, this heterogeneity is the consequence of a large increase of cooperativity in the 
vicinity of Pc>- 

4-3. Thermodynamics of supercooled water. 

Four scenarios have been proposed to explain the thermodynamics supercooled water. 
The stability limit scenario |62j hypothesizes that the limit of stability of superheated 
liquid water merges with the limit of stretched and supercooled water, giving rise to 
a single locus in the P — T plane, with positive slope at high T and negative slope at 
low T. The reentrant behavior of this locus would be consistent with the anomalies of 
water observed at higher T. As discussed by Debenedetti, thermodynamic inconsistency 
challenges this scenario [63] . 

The liquid-liquid critical point (LLCP) scenario [M] supposes a first order phase 
transition in the supercooled region between two metastable liquids at different densities: 
the low-density liquid (LDL) at low P and T, and the high-density (HDL) at high P 
and T. The phase transition line has a negative slope in the P — T plane and ends in a 
critical point. Numerical simulations for several models are consistent with this scenario 

[B1E3I66JE7]. 

The singularity-free scenario [39J focuses on the anticorrelation between entropy and 
volume as cause of the large increase of response functions at low T and hypothesizes no 
HB cooperativity. The scenario predicts lines of maxima in the P — T for the response 
functions, similar to those observed in the LLCP scenario, but shows no singularity for 
T > 0. 

The critical-point-free scenario [68] hypothesizes an order-disorder transition, with a 
possible weak discontinuity of density, that extends to P < and reaches the supercooled 
limit of stability of liquid water. This scenario would effectively predicts no critical point 
and a behavior for the limit of stability of liquid water as in the stability limit scenario. 

As showed by Stokely et al. [51], all these scenarios may be mapped into the 
space of parameter J and J a , of the model presented in the previous section, i.e. the 
coupling constants of the covalent component of the HBs and the coupling constant of 
the many-body component of the HB interaction, respectively. In particular, Stokely et 
al. showed by mean field calculations and numerical simulations that the absence of the 
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many-body component leads to singularity free scenario, while a large value of the many- 
body component with respect to the covalent component give rise to the critical-point 
free/stability limit scenario [33] . 

By using estimate of these parameters from experimental results, the authors 
predict a liquid-liquid phase transition at low temperature and high pressure ending in a 
second critical point C for water [51]. Therefore, following this prediction, the increasing 
fluctuations related to Cp, K T and «p of water under cooling are consequence of the 
liquid-liquid critical point C in the supercooled region of water. By approaching C, 
the correlation length £ of the HBs increases. In particular for any P < Pc, the critical 
pressure, there is a temperature TV where the correlation length £(P) is maximum. 
The locus TV(P), called Widom line, converges toward C with a negative slope in the 
P — T plane [HI EH]- By increasing P along the Widom line, £ increases and diverges 
at P = Pc- Therefore, the regions of cooperativity of HBs increase in size, leading to 
long cooperativity and, as a consequence, to larger heterogeneity in the dynamics as 
observed by Franzese and de los Santos [52] (discussed in the previous section). 

Before discussing in more details the features of these cooperative regions, is worth 
mentioning that recent theoretical and experimental results on water hydrating lysozyme 
proteins at very low hydration level (h = 0.3 g H^O/g protein) allow to explain the very 
low-T phase diagram of water monolayer, at T ~ 150 K at ambient pressure [70]. This 
investigation reveals that at low P two structural changes take place in the HB network 
of the hydration shell. One at about 250 K is due to the building up of the HB network 
[57], and another at about 180 K is consequence of the cooperative reorganization of 
the HBs. These two structural changes give rise to two dynamical crossovers in the HB 
correlation time and the corresponding experimental quantity, the proton relaxation 
time [70]. By increasing P, approaching Pc>, the two structural changes merge and 
at Pc 1 lead to diverging fluctuations associated to the liquid-liquid critical point, as 
discussed in a recent work by Mazza et al. |71j . 

5. Geometrical description of clusters of correlated HBs 

As discussed in the previous section, a water monolayer between hydrophobic plates 
separated by less than 1 nm, has a complex phase behavior at low T below the limit of 
stability of bulk liquid water. The same phase diagram compares well with experiments 
with water monolayer hydrating a complex substrate formed by proteins at low hydration 
level [SHI ED] • This can be understood if we admit that the main effect of the protein 
substrate at low hydration is to induce in the first layers of hydrating water a structure 
that is inconsistent with any possible crystal. Therefore, the substrate inhibits the 
crystallization, but does not inhibit the water-water HB formation. It is, therefore, 
interesting to understand how the region of correlated HB builds up and give rise to the 
cooperative rearrangement and the liquid-liquid phase transition. 

To this goal we follow an approach that has been validated during the last three 
decades to describe critical phase transitions. It consists in a percolation approach 
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elaborated in 1980 by Coniglio and Klein [72] for ferromagnetic systems and related to 
a mathematical mapping introduced by Fortuin and Kasteleyn [73J and to Swendsen- 
Wang [71] and Wolff [75] techniques for cluster MC methods. The Coniglio-Klein 
approach, called random- site- correlated-bond percolation, was extended ferromagnetic 
systems with many states [76] and spin-glass-like systems [77J [TSJ [7JJJ ED]- In particular, 
in Ref. [72] it was proved that the clusters defined following the rules of this specific 
type of percolation statistically coincide with the region of thermodynamically correlated 
variables. Moreover, in Ref. [78] it was proved that this result holds as long as the system 
has no frustration due to competing interactions. Since in the case considered here there 
are no competing interactions, we can follow the percolation approach to define clusters 
of water molecules with statistically correlated HBs. 

As described in Ref. [43J, we adopt the Wolff cluster MC algorithm [53J to study the 
cooperative regions and their length scale. Thanks to the fact that a cluster represents 
a region of water molecules with statistically correlated HBs, the algorithm allows to 
equilibrate the system at any T [53] . 

By definition, two variables and belong to the same cluster with probability 

p = mm{5 a ^ jV 1 - exp[-(J - Pv HB )/kT]} (22) 
if they belong to nearest neighbor molecules, or with probability 

p a = mm{6 aaaik , 1 - exp(-J CT /£;T)} (23) 

if they belong to the same water molecule i. 

The size of a cluster is given by the total number of variables belonging to the 
cluster. For each four in a cluster we have, on average, one water molecule in the 
cluster. The average linear size of finite (non percolating) cluster is, for the mapping 
discussed above, statistically equivalent to the correlation length of the HBs. Moreover, 
it is possible to prove [76] that each thermodynamic quantity, such as the compressibility, 
can be described in terms of an appropriate moment of the finite cluster distribution. 

By approaching the critical point C, we observe that the largest cluster percolates 
and its linear size becomes comparable to the system size. Under these conditions, the 
correlation length £ diverges. While away from C the distribution n(s) of cluster of 
linear size s decays as an exponential, n(s) has a power law decay near C (Fig. [2]) 
consistent with the theory [78J. 

From general considerations it is possible to show that n(s) ~ s~ T , where the 
exponent r is related to the fractal dimension Dp of the system r = 1 + d/Dp, and 
d = 2 is the embedding (euclidean) dimension. A preliminary estimate r ~ 2 suggests 
that the clusters of correlated HBs are compact with Dp ~ 2 |81] . Therefore, the 
mapping of the thermodynamic systems into a percolation problem allows us to give a 
geometrical description of the correlated HBs. 
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Figure 2. Distribution n(s) of clusters with finite size s formed by correlated 
hydrogen bonds (HBs) in a non-crystallizing water monolayer. Calculations are for 
P = 1.93 GPa ~ Pc'i the liquid-liquid critical pressure, and different values of 
temperature for system with N = 4 x 10 4 water molecules. For T — 173.84 K ~ Tq>, 
the liquid-liquid critical temperature, calculations (blue square points) decays as a 
power law n(s) ~ s~ T with r = 2.1 ± 0.1, as expected from theory near a critical point 
[75]. Consistent with the theory, we find that n(s) cannot be described by a power 
law decay away from the critical point. This is the case, for example, at T = 173.83 K 
(green circles) and T = 176.35 K (orange triangles). We find that at temperature 
far from the critical temperature, n(s) has an exponential decay, as expected. This is 
the case, for example, at T = 833.08 K (light blue square) and T = 163.80 K (black 
triangles). Continuous line is the power law fit, while dashed line are exponential fits. 



6. Free energy landscape analysis 

The percolation approach allows us to adopt a cluster MC dynamics that is very efficient 
at low T [53]. Therefore, we can equilibrate the system at low P around the T W (P), 
the temperature of maximum correlation length £, and around the temperature Tll 
of liquid-liquid coexistence at high P, and calculate the free energy landscape for the 
system. 

By definition, the Gibbs free energy is 

G/k B T=-\n^>(H,p), (24) 

where £P(H,p) is the density of states with enthalpy H G [H,H + 5H] and density 
p E [p, p+Sp], with 8 H and Sp infinitesimal increments. In Fig. [3] we show G as a function 
of energy per particle E/N and the density p with two equivalent minima straddling 
the line of phase transition. The two minima, equivalent within the numerical precision, 
correspond to two coexisting phases with different density and energy. The one at higher 
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Figure 3. Gibbs free energy for a water monolayer with N = 4 x 10 4 water molecules 
at P = 1.98 GPa and T = 158 K. The two minima, one at low energy and density, and 
the other at high energy and density, represent respectively the LDL phase and HDL 
phase of the system. The projection of G on the plane E/N — p shows that there is 
a linear relation between the accessible energies and density for the coexisting states. 
The same value of the minima marks the coexistence of the two phases. 

density and energy represents the HDL phase. The other at lower density and energy 
represents the LDL phase. Approaching C the two minima get closer and the density 
separation disappears, as expected at the critical point. These results are consistent 
with the mean field free energy analysis of Ref . [53l 182] , where the Gibbs free energy is 
calculated as function of the HB order parameters relevant at the structural transition 
at high P and low T. In the mean field analysis the minima of G are separated at high 
P, but merge for P approaching Pqi. All these results are consistent with the behavior 
of Cp, Kt and ap whose maxima move to lower T as P is increased [561 SHI EE1 158] . 
The loci of the maxima of Cp, and ap merge in the vicinity of C and the amplitude 
of their maxima increases approaching C . 

7. A model for protein in water 

In the previous section we define a coarse-grained model for a water monolayer and we 
show that the model compares well with experiments for protein shell-water and that it 
predicts a complex phase diagram for low T, below the limit of stability of bulk liquid 
water, and high P. As discussed in the introduction, under these conditions folded 
proteins are destabilized. Following our discussion about how could be relevant to take 
into account the HB free energy to explain the lost of stability of folded proteins, It is 
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intriguing to test if the proposed water model could give insight into the mechanism of 
unfolding. 

To this goal we modify the water model to introduce the effect of the protein-water 
interface. For sake of the simplicity, we will limit our discussion to the case of a single 
protein embedded into a water monolayer. Although far from the complex studies of a 
protein embedded into bulk water, the model gives instructive results. 

The simplest protein that we can consider is a hydrophobic homopolymer, 
schematized as a self avoiding chain (Fig. Q. Following the discussion by Muller |15j . we 
require that, consistent with experiment, water molecules in contact with a hydrophobic 
monomer have larger decrease of enthalpy upon HB formation. Also consistent with 
Muller-Lee-Graziano discussion [25], the fraction of broken HBs at the hydrophobic 
interface is larger than the fraction of broken HBs in the bulk. 

The first requirement is achieved by adding a term to the water Hamiltonian 



Eq. (18) 



= -A J ^ "'"/VAv- (25) 

<ij>s 

where the sum is taken over nearest neighbor water molecules in the protein hydration 
shell (Fig. [5]), and A > is an adjustable parameter accounting for the larger enthalpy 
decrease for HBs in the hydration shell. Hence, for HB formed between water molecules 
in the shell, the enthalpy variation is — J(l + A) + Pvhb, and the total enthalpy for 
protein into water is 

H tot = H + J%, (26) 



where H is given by Eq. (18). 

The second requirement of Muller-Lee-Graziano approach, i.e. a larger number 
of broken HBs at the interface, is achieved by volume exclusion. Once a cell of our 
system is occupied by a protein monomer, it cannot be occupied by a water molecule. 
Therefore, a water molecule, in the hydration shell cannot form the HB in the direction 
of the monomer and looses at least one HB (it can loose more if it has more monomers 
as nearest neighbors, as shown in Fig. [1]). 

In the following section we will define the algorithm adopted to generate protein 
equilibrium configurations. To this propose we follow a MC procedure that mimics the 
dynamics at large time scales. 



7.1. Monte Carlo algorithm 

We perform MC simulations in the NPT ensemble. In every MC step we choose a 
cell at random. If it is occupied by a water molecule, we change randomly one of its a 
variables. If it is occupied by a monomer and if the monomer is in a corner configuration 
(Fig. [5]^.) then we swap its position with the position of the water molecules in the cell 
in the opposite corner. By doing this, we keep the inter-monomer distances constant. 

If the cell, picked at random, is occupied by a monomer not in a corner 
configuration, no displacement is performed because it would change the inter-monomer 
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Figure 4. Example of configuration of a homopolymer in the coarse-grained model. 
Each cell is occupied either by a water molecule (white and gray cells) or a hydrophobic 
homopolymer monomer (cells with a full black circle). The gray cells represents the 
sites occupied by shell water. The enthalpy gain for HB formation between shell water 
molecules is larger than that between bulk water molecules, according to the Eq. 25 
Shell water molecules cannot form hydrogen bonds with nearest neighbors hydrophobic 
monomers. 



distance. This limitation is introduced in order to avoid in the free energy any term 
accounting for the elastic energy of the homopolymer. The effects of this elastic 
contribution are for the moment outside of the scope of the present work. 

Finally, as in the cases discussed in the previous section, to keep the pressure of 
the system constant, every iV random changes of the cell variables (where N is the 
total number of cells in the system), we attempt to rescale all the system volume by a 
factor that is tuned in a way to guaranty 50% of acceptance ratio. All the MC moves 
described above are accepted or rejected according to the Boltzmann factor associated 
to the enthalpy change caused by the move. 

In order to study the folding-unfolding process of the proteins we calculate the 
number of contact points iV C p ts as illustrated in Fig. (|5]). In this calculation we do not 
count the monomers that are adjacent along the homopolymer. 

8. Preliminary results for hydrated homopolymers 

We study a system with N water molecules and a hydrophobic homopolymer chain 
with iV m monomers. In our preliminary simulations we used N = 650 or N = 1000 
and iV m = 12 or N m = 50. The parameters are chosen, for consistency, as in previous 
analysis [S3]: e = 5.8 kJ/mol, J = 2.9 kJ/mol, J a = 0.29 kJ/mol, v = hr 2 Q) h = 7 A, 
v-rb/vq = 0.5 and q — 6. We choose A = 0.7 for the larger decrease of enthalpy at the 
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Figure 5. A) Monomer 2 is in a corner configuration and can be displaced from 
the configuration on the left to the configuration on the right and vice versa. B) 
Homopolymer configuration with two contact points, indicated with dotted lines. C) 
Homopolymer configuration without contact points. 

hydrophobic interface. 

Our results display a non monotonic behavior of -/V cpts as function of T at low P. 
At high pressure we observe an region in the P — T plane where the number of contact 
points is at least 51% of the maximum possible number. By definition, we consider 
these configurations as belonging to the set of folded states. We observe that the region 
of folded states is included within a larger region in the P — T plane where the number 
of contact points is at least 49% of the maximum possible number. By definition, we 
consider those configurations as members of the set of state representing the molten 
globule [831. 

We find that the region of folded states has an elliptic shape that resembles the 
theoretical prediction (Fig. [T]). In particular, we observe that a folded protein unfolds 
upon cooling, giving rise to the cold denaturation process. It also unfold by increasing 
the pressure as expected by pressure denaturation (Fig. [6J. Since our stability region 
is at high P, we are also able to observe the unfolding by decreasing the pressure, a 
phenomena that is predicted by general theoretical considerations, as discussed in the 
introduction. We also find that the axes of the elliptical stability region are tilted as 
expected (Fig. [I]). 

9. Conclusions 

The behavior of supercooled water is still under debate and the presence of a second 
critical point C could be relevant to understand how the structure of liquid water 
changes around proteins and how affects protein properties. Experiments of water 
confined in nano-structures offer the possibility to access a range of temperatures where 
bulk liquid water would not be stable and would form ice. Hence, confinement allows 
to study water under conditions that are relevant in biological systems. 

Despite the growing interest of the scientific community in water at hydrophobic 
and hydrophilic interfaces, it is still unclear how the interaction with the confining 
surfaces affects the thermodynamics of water. For example, recently Strekalova et al. 
[EH] observed that the fluctuations of supercooled water confined into a hydrophobic 
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Figure 6. Typical configurations of folding-unfolding of a coarse-grained protein 
suspended in water at different temperatures T and high pressures P. The protein is 
represented as a fully hydrophobic chain (in white), surrounded by water molecules 
(turquoise background). We use different color for water to indicate the different 
orientations of the HBs. (a) At high P and high T, the protein unfolds and the number 
of HBs (colored sticks) of surrounding water is small, (b) At the same pressure but 
lower T, the protein collapses in a molten globule state, (c) At lower T the protein 
folds, while the surrounded water has a large number of HBs. (d) At much lower T we 
observe cold denaturation of the protein when the number of HBs is largely reduced 
(zero in the configuration represented here), (e) At higher P the denaturation occurs 
at higher T, and the mechanism of unfolding seems to be dominated by the reduction 
of HBs also under these conditions. 



porous material are drastically smaller than those of bulk water. They found that the 
response functions Cp, ap and Kt are largely reduced as a consequence of the interaction 
with the porous medium. An extreme consequence of this change is the disappearing 
of the liquid-liquid phase transition at high pressures [83]. Therefore, further work is 
necessary to clarify the many issues related to the dynamics and thermodynamics of 
water at the interfaces. 

Here we presented a coarse-grained model for a monolayer of water and its extension 
to the case of solvated proteins. The model takes into account the cooperativity between 
HBs and has been studied by simulations and mean field calculations. Previous results 
about the phase diagram, the diffusivity properties, the response functions Cp, ap and 
Kt of the model and the connection of these quantities with the HBs dynamics are in 
agreement with experimental results and validate the model. 

We adopted this model in the context of protein folding. For the sake of 
simplicity we consider the case of a protein schematized as a self-avoiding hydrophobic 
homopolymer. Following Muller analysis [TH], we assume that the network of HBs is 
perturbed by the presence of hydrophobic solute with large size. 

Our preliminary results reproduce hot, cold and pressure denaturation as well as 
the existence of intermediate states (molten globule). We find that the stability region 
for folded protein has the theoretically expected elliptic shape in the P — T plane. Fur- 
ther work is in progress to elucidate the relevant mechanism ruling protein stability. 
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